Slow energy relaxation of macromolecules and nano-clusters in solution 
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Many systems in the realm of nanophysics from both the living and inorganic world display slow 
relaxation kinetics of energy fluctuations. In this paper we propose a general explanation for such 
phenomenon, based on the effects of interactions with the solvent. Within a simple harmonic model 
of the system fluctuations, we demonstrate that the inhomogeneity of coupling to the solvent of 
the bulk and surface atoms suffices to generate a complex spectrum of decay rates. We show for 
Myoglobin and for a metal nano-cluster that the result is a complex, non-exponential relaxation 
dynamics. 
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Matter from both the living and the inorganic worlds 
displays unusual structural and dynamical properties 
when it is reduced to small objects of nano-metric scale 
. One interesting feature that has recently been arous- 
ing interest is the relaxation to thermal equilibrium of 
local or distributed temperature fluctuations. Heat dissi- 
pation experiments in metallic nano-clusters in solution 
have revealed that the approach to equilibrium follows a 
strctched-exponential law 0- Similarly, experiments on 
proteins have shown that their relaxation after excitation 
may also be highly non-exponential 0, 3 ■ 

Slow energy relaxation is naturally associated with the 
presence of a hierarchy of relaxation times. Neverthe- 
less, the origin of such hierarchy is not always easy to 
pinpoint. A common assumption is that multiple time 
scales arise because of the roughness of the energy land- 
scape, i.e. lots of local minima of different depths whose 
sequence en route to equilibrium determines a cascade 
of decay times. However, while such a picture would 
be compatible with our current understanding of protein 
folding dynamics j^, such mechanism is unlikely to ap- 
ply to metallic nano-clusters, which are regular stacks of 
identical atoms, and whose energy landscape close to the 
ground state is likely to be mostlyfeatureless. Yet, the 
relaxation dynamics in solution of proteins and nano- 
clusters share enough similarities to raise the question 
whether a common mechanism could be at work in a 
broad class of nano-systems. 

The physiological environment of proteins is a viscous 
aqueous solution, while metallic nano-clusters are be- 
lieved to be promising markers to study the behavior 
of cells and living tissues, so that an aqueous solution is 
also the environment in many of their applications. In 
this paper we propose that a natural, simple and uni- 
fying explanation for the spectrum of decay times for 
nano-particles in solution can be found in the different 
coupling to the solvent of surface and bulk atoms. 

The idea is that energy is dissipated to the environment 
only by surface atoms, whereas bulk atoms only exchange 



energy with each other with very little or no dissipation. 
To investigate the effects of the damping inhomogeneity, 
we model the systems under investigation as elastic net- 
works. The interatomic interactions are simple harmonic 
potentials V{n,rj) = /2{\rt - rj \ - {r^ - r^ol)"^ , where 
ri is the position of atom i, fio hs equilibrium position 
and kij is the interaction stiffness between atoms i and j. 
In the case of metallic nano-clusters this modelization is 
the usual, textbook description, which correctly predicts 
most of the general features of crystal vibrations [q] . The 
application of elastic network models to proteins is more 
recent 0, since it had been customary to assume that 
proteins are characterized by complex energy landscapes. 
Yet, it has been realized that most features of the large- 
and medium-scale dynamics of proteins close to their na- 
tive state can be successfully reproduced by simple har- 
monic interactions between amino-acids H H [13. In 
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Figure 1: Structure of Myoglobin (pdb code 1A6M) and of 
a model icosahedral metal cluster with N = 309 atoms. The 
protein structure is color-coded according to the fraction of 
surface accessible to the solvent at each site - from blue (low 
accessibility) to red (highest exposure). 
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these coarse-grained models represents the position of 
the a-carbon of the i-th amino-acid, r^o its position in the 
native state as determined from X-ray crystaUography or 
Nuclear Magnetic Resonance, and kij can take different 



= k e{\fio - fjo\ - re) or 
where Vn is a suitable cut- 



functional forms, such as fcy 
hj = k cxp(-|r',o - rjolV^c) 
off (or typical) interaction distance that tunes the overall 
connectivity of the structure, and k a phenomenological 
strength constant. 

Treating the systems as elastic networks is indeed the 
simplest hypothesis. Yet, recent rigorous findings about 
energy relaxation in one and two dimensions suggest that 
even linear distributed systems display highly non-trivial 
relaxation properties pH IT^ . We thus start up illus- 
trating the solution to the following problem: a chain of 
N equal masses connected by springs of equal strength, 
coupled to a heat bath at its edges. The chain has free 
ends and at t = is at equilibrium at a temperature T,,. 
Modeling the system a la Langevin, only the first and 
last particles are damped and subject to randomly fluc- 
tuating forces due to their interaction with the solvent, 
according to the equations 



mxi = k[xi-i — 2xi + Xi-i-i) i = ^ 
mx'i = k{x2 — xi) ~ jxi + T]i 
mxN = k{xN-i ~ xn) - jx'n + Vn 



,N-1 



(1) 



where 7 is the viscous friction coefficient and r]k{t) (fc = 
l,N) is a Gaussian, delta-correlated white noise, whose 
standard deviation is fixed by the fluctuation dissipation 
theorem. We want to study the system relaxation to a 
temperature T < Tq. 

The time behavior of the system only depends on the 
temperature difference AT = T„ — T because of the ab- 
sence of local energy barriers in the energy landscape. 
This property shall be proved rigorously later on in the 
framework of the Fokker-Planck (FP) approach. As a 
consequence, we may solve the problem as that of a de- 
terministic dissipative system of energy NkgAT damped 
at its edges. Such problem is amenable to a perturbative 
treatment in the limits of small and large damping, re- 
sul ting in a decay spectrum of the linear modes given 
by 
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where 7 and 0o are the damping strength and band- 
edge frequency, respectively, and / = 2/N is the sur- 
face fraction. The superposition of decay constants of 
the form Q allows to calculate the energy relaxation an- 
alytically, revealing a crossover in the system from the 
initial exponential decay exp(— 7ot) (single fastest relax- 
ing mode) to a power-law relaxation of the type 
(mte(7raied regime) [T]] |. Accordingly, the solution to the 



original problem may be written as 

{E{t)) - {E{^)) 



{E{Q)) - {E{oo)) 



I .9(0)e-2^W*d0 = Xo(7oi) 




(3) 



where g{i-o) is the density of modes and Iq is the 
exponentially-modified zero-order Bessel function. The 
same calculation may be repeated in higher dimensions, 
yielding 
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where d is the spatial dimension, the parameter 70 still 
carrying the information on the linear surface fraction 
f/d. In the rest of the Letter we show that a simi- 
lar crossover from an exponential to an integrated, non- 
exponential regime also occurs for more complex geome- 
tries and connectivities. In order to treat the general 
case, we numerically compute the exact solution of the 
FP formulation of the problem. 

In the harmonic approximation, the total potential en- 
ergy of a system of size N can be written as a quadratic 
form 
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{X ~ Xy K {X - X") 



(5) 



where X = {ri,^, ri,j^, ri,^, . . . , tat^^,, tat^^, rjv,z}, X° the 
same vector at equilibrium and the "contact" matrix K 
is simply the Hessian of the potential energy function 
evaluated at the equilibrium structure (for the sake of 
simplicity, we have set all masses equal to one) . Accord- 
ingly, the FP equation for the probability distribution in 
phase space takes the form [l^ 
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P{Y,t\Ym (6) 



where Y ~ {X ~ X° ^ X) is the 6A^-dimensional vector of 
displacements and velocities, and the matrices A and B 
are given by 




^ksT 




(7) 



The information on the coupling to the solvent of indi- 
vidual atoms is contained in the diagonal matrix F. The 
latter has the form F^ = jSijSi, where the vector S fixes 
the fraction of surface exposed to the solvent by each 
particle (0 < 5^ < 1, i = 1, 2, . . . , TV). The parameter 7 
specifies the overall strength of the viscous force. 
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Figure 2: Relaxation to equilibrium in Myoglobin (a) and in an icosahedral cluster (b). Symbols are plots of T^{t) (see text) 
as calculated through Eq. 1131 . Squares: under-damped regime (7/m = 0.27 ps"^ (a) and 7/m = 4.27 ps"^ (b)). Circles: 
over-damped regime (7 = 27 ps^^ (a) 7 = 427 ps~^ (b))- The dashed lines are plots of the straight lines ©(t) — 2fcS'yt. The 
dotted and dot-dashed lines are stretched exponential laws with exponent a ~ 0.22 and a ~ 0.24, respectively (a) and a ~ 0.35 
and a ~ 0.48, respectively (b). The solid lines are plots of formula with 70 = (/cfi/dcff)7 and effective dimension dcff ~ 1.45 
(a), and d^B « 0.83 (b). In the insets we plot {E{t)) — {E{oo)) in units of ksT (symbols). The solid lines are the same as in 
the main panels. Parameters are To ^2T, N = 194, f^s = 0.242, = 8 A (a) and iV = 309, /off = 0.217, = 6 A (b). 



The solution of Eq. © is the multivariate Gaussian 
distribution Q 

p(y,t|y(o)) = {2TT)-^^\dctc{t)\-^/^x 

exp - G{t)Y{0)fC-\t)[Y - G(i)r(0)]| (8) 

{Yi{t)Yj (t)) is the block correlation ma- 
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and G(t) is the propagator matrix. The latter can 
be evaluated at any time in terms of a normalized bi- 
orthogonal set of left and right eigenvectors of the matrix 
A (sometimes referred to as Langevin modes) as 

G{t) = ^fle^'^-i (10) 

where A is the diagonal matrix of the eigenvalues of A, 
and ^'ij and ^'l are the matrices of right and left eigen- 
vectors, respectively. It is easy to show that the evolution 
law for the correlations reads, in matrix form, 

C(t) = C{oo) + G{t)[C{Q) - C{oo)]G^{t) (11) 

where we have introduced the equilibrium correlation ma- 
trix 



C(oo) = ksT 







(12) 



K^^ being the generalized inverse of matrix K. 

In order to calculate the energy relaxation to equilib- 
rium, we first diagonalize matrix A and then evaluate 
the correlation matrix Hll|l as a function of time with the 
aid of Eq. H1U|) . The expression for the total energy then 
finally reads 
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{E{t))^-TT[G^^{t)+KCxx{t)] 



(13) 



To illustrate our analysis, let us examine the relaxation 
to equilibrium of a typical globular protein. Myoglobin, 
and of a model Au nano-cluster with icosahedral symme- 
try ^3 (see Fig. We estimate the eflFective fraction 
of surface exposed to the solvent at each site (the vector 
S) through standard solvent-accessible surface areas. As 
initial conditions, we set the specific potential energy at 
its equilibrium value 3/2kBT, while each particle is given 
the kinetic energy "i/iksTo {T„ > T). This amounts 
to taking in Eq. ^ Cxx{0) = Cxx{oo) = k^TR'-^ 
C^xiO) ^ C'^x(oo) = 0, and C_^^(0) = ksTo\^. Inci- 
dentally, the latter assignments, together with Eq. Hll|l . 
also prove that the relaxation process only depends on 
the diflFerence AT = Tq - T. 

Our results are summarized in Figs. |21 (a) and (b), 
where we plot the quantity Vlt) defined as 



P(i) = -log 



{E(t)) {E{^)) 



{Em {E{^)) 



(14) 



This representation is most useful in identifying 
crossovers from an exponential law (straight line) to a 
slower decay regime, such as a stretched exponential of 
the form exp[— (t/r)'^] (power law with cr < 1). 
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Remarkably, the two systems under scrutiny display 
the same complex behavior of the energy decay. Ini- 
tially, the relaxation is exponential, which is the first 
signature of the linear damping on surface atoms. Ac- 
cordingly, Vit) = 2fcsjt/m, where fcs = Y^t^i/N is 
the total effective surface fraction of the system. The 
first crossover occurs in correspondence with the onset of 
energy transfer from kinetic to potential (arrows marked 
Ci in the plots), which occurs on a time scale of the order 
of the inverse of the maximum frequency. This crossover 
is independent of the damping rate 7. At this stage the 
potential energies of all modes have been transferred a 
fraction of the initial kinetic energy, and thus the relax- 
ation process can be described in terms of Eq. |0J. In 
fact, it turns out that the decay curves in the under- 
damped regime can be extremely well approximated by 
formula with 7,, = (/off/rfcff)7, where the linear sur- 
face fraction is naturally obtained by allowing for an ef- 
fective dimension deft (see insets in Figs.|21(a) and (b)). 
One-parameter fits of the exact solutions yield deff ~ 1.45 
and defi ~ 0.83 for Myoglobin and the nano-cluster, re- 
spectively, which agree with the sizeable surface fraction 
of the systems. 

In accordance with the prescription of formula , at 
a time of the order des/ fesil/m) the systems undergo a 
second crossover to an integrated regime, which reflects 
the superposition of the time scales of all Langevin modes 
(arrows marked C2 in the figures). This second crossover 
is by all means captured by the transition from expo- 
nential to a power law implicit in formula (0J. In the 
over-damped case, the fast damping of the surface atoms 
increases the time scale for the dissipation of the bulk 
energy. As a consequence, the first stage of surface relax- 
ation stretches until the system reaches the crossover to 
the integrated regime, thus merging the two crossovers. 

In our chosen time span, the long-term behavior of the 
decay curves is best approximated by a stretched expo- 
nential (see again Figs. 12)1 . To this concern, it is however 
important to stress that asymptotically the decay laws 
shall eventually bend again to a pure exponential. This 
is a finite-size effect that reflects the progressive depletion 



of modes, until only the slowest one is left with signifi- 
cant energy. Consequently, its decay constant j„ sets 
the time scale for the last exponential stage exp(— 27jvt), 
which in both systems is of the order of 10^ in natural 
time units (00 = 2). This effect may make it hard to 
distinguish between a stretched exponential and a slow 
transition from power-law to exponential. 

In this Letter we have shown that the inhomogeneity of 
coupling to the solvent of the surface and core atoms nat- 
urally produces a hierarchy of relaxation times, and hence 
complex energy relaxation, in systems as un-suspect as 
harmonic bead-spring networks. It can be argued that 
non-linear effects could also affect the nature of relax- 
ation dynamics [l6[ . However, within our treatment, at 
least for small polynomial non-linearities, an-harmonic 
effects are not expected to introduce dramatic modifica- 
tions. In fact, the fast damping of low-frequency modes 
with respect to band-edge ones rapidly hinders effective 
inter-mode energy exchange. This phenomenon has been 
well documented in one- and two-dimensional lattices 
with quartic anharmonicity [Til Il2l | , where the crossover 
from exponential to a collective, power-law regime has 
been observed during relaxation irrespective of the an- 
harmonicity strength. It is however interesting to recall 
that relaxation in non-linear systems is under certain cir- 
cumstances associated with spontaneous energy localiza- 
tion 0, This phenomenon is definitely worth at- 
tention for example in the context of protein dynamics, 
as it may play a role in relaxation and redistribution of 
local energy fluctuations, such as after hydrolysis of ATP 
molecules. 

We stress that our results are qualitatively consistent 
with a number of experimental observations 0, 0| • How- 
ever, further experimental work would be necessary in 
order to discriminate between our model and alternative 
explanations, such as ruggedness of the system energy 
landscape Q. 
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